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Weyl semimetals are gapless quasi-topological materials with a set of isolated nodal points form¬ 
ing their Fermi surface. They manifest their quasi-topological character in a series of topological 
electromagnetic responses including the anomalous Hall effect. Here we study the effect of disor¬ 
der on Weyl semimetals while monitoring both their nodal/semi-metallic and topological properties 
through computations of the localization length and the Hall conductivity. We examine three dif¬ 
ferent lattice tight-binding models which realize the Weyl semimetal in part of their phase diagram 
and look for universal features that are common to all of the models, and interesting distinguish¬ 
ing features of each model. We present detailed phase diagrams of these models for large system 
sizes and we find that weak disorder preserves the nodal points up to the diffusive limit, but does 
affect the Hall conductivity. We show that the trend of the Hall conductivity is consistent with an 
effective picture in which disorder causes the Weyl nodes move within the Brillouin zone along a 
specific direction that depends deterministically on the properties of the model and the neighboring 
phases to the Weyl semimetal phase. We also uncover an unusual (non-quantized) anomalous Hall 
insulator phase which can only exist in the presence of disorder. 


I. INTRODUCTION 

The topological classification of symmetry protected 
fermionic systems has added a new level of complexity 
to the study of quantum phases of matter. The clas¬ 
sification scheme was originally considered for gapped 
phases^^^, motivated by the discovery of symmetry pro¬ 
tected topological phases with discrete anti-unitary sym¬ 
metries^’®. Over the past few years there has been a 
tremendous interest in understanding various topologi¬ 
cal phases in gapless materials as well, both experimen¬ 
tally®^^^ and theoretically®®^^® . Aside from graphene 
in two dimensions, the Weyl semimetal (WSM) is the 
most celebrated example of a quasi-topological gapless 
phase in three dimensions. Its band structure is gapped 
everywhere in the Brillouin zone (BZ) other than few iso¬ 
lated points, i.e., Weyl nodes, where the conduction and 
valence bands touch to form a doubly-degenerate point. 
The Weyl nodes are locally stable, topological objects in 
momentum space®^, and they each carry a chirality quan¬ 
tum number (Berry phase magnetic charge). Owing to 
the presence of these nodes, a WSM can display exotic 
properties such as surface Fermi arcs®^®®’®®’®®’'®®, negative 
magnetoresistance®’®®’®®, a magnetic torque anomaly®®, 
and an anomalous Hall effect®’®®’®®^'®®. 

Weyl nodes are believed to be robust against pertur¬ 
bations unless discrete translation symmetry or charge 
conservation symmetry is broken, both of which would 
allow the nodes to mix with each other and open a gap. 
In this paper, we investigate the robustness of WSM 
phases after breaking translational symmetry via real- 
space quenched disorder. Stability of a single Weyl node 
against disorder has been a subject of great effort using 
continuum limit calculations®®’®®^®®, as well as numerical 
simulations®®^®®. The consensus is that a Weyl node (be¬ 
sides non-perturbative rare region effects®®’®®) remains a 
node, with vanishing density of states, up to some finite 
critical disorder, at which it transforms into a diffusive 


metal with a finite density of states at the nodal point. 

Unlike these previous studies, the focus of our arti¬ 
cle is on lattice tight-binding (TB) models with pairs of 
Weyl nodes. Considering TB models is motivated by the 
following reasons. First, a single Weyl node cannot ap¬ 
pear in isolation in a solid state material according to 
the Nielsen-Ninomiya no-go theorem®^. Second, some of 
the gapped regions in the WSM phase can be thought 
of as layered topological phases (stacks of Chern insu¬ 
lators) in their own right. The gapped regions lead to 
their own experimental consequences (anomalous Hall ef¬ 
fect, surface states, etc.) and may impact the outcome. 
Hence, other observable properties such as Fermi arcs or 
the anomalous Hall effect (AHE) need to be taken into 
account in order to determine the nature of a disordered 
WSM phase. Measuring topological responses requires 
the full lattice model to include all non-local (in momen¬ 
tum space) contributions and cannot be done using the 
low energy effective theory alone. Third, it is more nat¬ 
ural to address the stability of a WSM with reference 
to its neighboring phases in the full phase diagram of a 
given lattice model. This is a clear way to identify which 
phase(s) the WSM transitions to as the system becomes 
more disordered. 

In order to quantify our investigation of the stabil¬ 
ity of WSM phases we need to keep track of quantities 
which are special to the WSM and thoroughly charac¬ 
terize it. In other words, one must check simultaneously 
what disorder does to the topological character of the 
occupied states as well as the semi-metallic behavior of 
the Weyl nodes. To monitor the topological properties, 
we study the AHE by computing the Hall conductivity 
Oxy, and to establish the semi-metallic behavior we cal¬ 
culate the localization length. We find that the usual 
expectation that the disorder will mix the Weyl nodes 
and tend to destroy the WSM phase is not always the 
case. Indeed there are instances where disorder makes 
the WSM phase stronger. There are several interesting 
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recent works®®^®® that have also studied TB models of 
Weyl/Dirac semimetals. In these works they have been 
interested in critical scaling properties near the metal- 
insulator transition point evaluated by computing den¬ 
sity of states, localization length, and/or the longitudinal 
conductivity. Our work is markedly different from these 
studies in the sense that we consider both semi-metallic 
and topological aspects of disordered WSMs. A recent, 
independent study^® in parallel to our work has also pur¬ 
sued this direction, and has mapped out a phase dia¬ 
gram for a disordered multi-layer Chern insulator model, 
albeit with a different type of disorder. In that work 
the phase boundaries are also found to match with the 
self-consistent Born approximation (SCBA) calculations 
in the weak disorder limit. In our work, we consider 
three models, the first of which is similar to the one par¬ 
tially studied in Ref. 70 (where the results overlap they 
are in complete agreement), and we attempt to recog¬ 
nize the common features among them. In explaining 
our observations, we try to adopt a view based on the 
generic properties of models (in terms of the distribu¬ 
tion of Chern numbers in the BZ) rather than trying to 
solve the SCBA explicitly for each model. Moreover, in 
terms of methodology, we compute Uxy using the much 
more efficient Green function approach^^^^® compared to 
the exact diagonalization method used in Ref. 70. The 
Green function method lets us carry out our calculations 
for large system sizes where exact diagonalization is un¬ 
wieldy. 

As mentioned, one of the essential goals of our work 
is to identify which properties and trends of disordered 
WSM lattice models are universal. For this purpose, 
we study three different lattice realizations of WSMs. 
These TB models contain at least one WSM phase in 
their phase diagrams: (i) a multi-layer Chern insula¬ 
tor (Cl) where layers have Chern number C = ±1, (ii) 
a double WSM model with layers with Chern number 
C = 2, and (iii) a 3D, four-band Dirac model tuned 
near the critical point between the topological insulator 
(TI) and trivial insulator phases, subject to a Zeeman 
splitting, (i) and (iii) support Weyl nodes with Berry- 
monopole charge ±1 while (ii) is an example of double- 
WSM having Weyl nodes with Berry-monopole charge 
±2. In summary, our results show that the Weyl nodes, 
regardless of their charge and model Hamiltonian, survive 
at weak disorder and effectively move either towards or 
away from each other. We find that the movement di¬ 
rection is determined by the sign of a quadratic term 
in the in-plane momenta normal to the line separating 
the Weyl nodes. Moreover, we find that a normal in¬ 
sulator (NI) may be transformed into a WSM at finite 
disorder strength. In this case Weyl nodes nucleate at a 
finite disorder strength and give rise to a finite AHE re¬ 
sponse. Inversely, disorder can hybridize the Weyl nodes 
and hence transform a WSM to a 3D Cl (or TI) with 
a bulk (or surface respectively) AHE. We show that all 
these observations are consistent with the SCBA analy¬ 
sis of disordered CIs^^. Overall, the effect of the SCBA 


is to shift the phase boundaries in each phase diagram 
as a function of disorder strength. In multi-layer CIs, 
models (i) and (ii), the WSM region of the phase dia¬ 
gram is maintained up to the diffusive limit; whereas in 
model (iii) a novel disorder-driven insulating phase with 
non-zero bulk AHE emerges at the top of WSM region 
before merging into the diffusive metal. In all cases, the 
AHE is continuous as the WSM region transitions into 
the diffusive metal (DM) and smoothly decreases to zero 
on the diffusive side as we reach the Anderson localized 
regime. This decreasing trend in the AHE is a direct 
consequence of diffusive extended bulk states which al¬ 
low scattering between counter-propagating edge modes 
on opposite sides through the bulk. 

In Sec. II we start by briefly introducing our methods, 
next we discuss the phase diagrams of the above models 
in three subsequent Secs. II A, II B, and II C, respectively; 
finally, we conclude our article in Sec. HI. 


II. TIGHT-BINDING MODELS FOR THE 
WEYL SEMIMETAL 

In this section, we introduce three tight-binding mod¬ 
els and study their properties as disorder is added to the 
system. Throughout this paper, the disorder is treated 
as a random on-site energy: 

^dis — ^ ^ ; (^■^) 

X 

where Cx and c^. are electron destruction and creation op¬ 
erators (and may have extra spin/orbital indices) and 
is a random number uniformly distributed in the range 
[—IT/2, IT/2] where IT is the disorder strength. For each 
model we consider we will present a phase diagram where 
each phase is characterized by the Hall conductivity 
and the localization length. In order to calculate the lo¬ 
calization length, we solve the standard recursive equa¬ 
tion by viewing the system as a stack of layers along the 
z direction^®^^^. The localization length is determined 
by ^ = I/k in terms of the smallest positive Lyapunov 
exponent k for a quasi-one-dimensional system of cross 
section x Ly = L x L and length Lz using the resol¬ 
vent (Green) function method which does not require any 
additional condition on the invertibility of the inter-layer 
hopping matrix. In our calculations, we take disorder en¬ 
semble averages for system sizes of L = 8,10,12,14 and 
Lz = 10^. We should note that the systems are quasi- 
one dimensional {Lz 3> L) which means that they are al¬ 
ways localized or equivalently the transfer matrix always 
decays exponentially with L^, i.e. T oc ex-p{—2Lz/^). 
However, one can still distinguish metallic phases from 
insulating phases by looking at how ^ behaves as the 
cross section is made larger. Let the normalized localiza¬ 
tion length be A = ^/L. Therefore, in a (semi-) metallic 
phase A is proportional to L which implies ^ oc (or 
the number of conducting channels), whereas in an in- 
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sulating phase A is decreasing as 1/L which means ^ is 
almost constant. Any critical point is characterized by a 
scale invariant A where all curves near this point fall onto 
a single-parameter scaling ansatz. We use these facts to 
determine the boundaries between various phases in our 
phase diagrams. We impose the periodic boundary con¬ 
ditions across each slice along x and y directions. 

For a clean WSM with two nodes at ±fco along the 
^-direction, the Hall conductivity in the normal plane is 
known to be 


^ xy — 


1 




Tch 


( 2 . 2 ) 


where a‘^y{kz) is the Hall conductance of a 2D layer at 
kz in the 3D BZ^®. In the rest of this paper, we set 
e = h = To compute a^y in a disordered media, we 
use the Green function method developed in Refs. TI¬ 
TS. The Green function approach is very efficient for 3D 
calculations since in this method the 3D sample is di¬ 
vided into 2D slices normal to a given direction (e.g., the 
x-direction) and the response function is computed itera¬ 
tively as more layers are added through a set of recursion 
relations. The idea is to write the Hall conductivity as 


~ T T T 


Tr 








(2.3) 


for a 3D sample with slices each containing LyL^ 
sites. Here is the retarded (advanced) Green’s func¬ 
tion connecting the i-th and j-th layers, and the position 
operators are denoted by x and y. Gfj is computed at 
Z = y, + ij where y is the chemical potential and the 
phenomenological lifetime parameter 7 is introduced to 
avoid singular behavior. Ultimately the limit 7 —>■ 0 is 
assumed. is evaluated at the (n -I- l)-th step 

(for n + 1 layers) from the Green function GfG'> (for n 
layers) by updating the self-energy with the contribution 
from the {n + l)-th layer. The efficiency of this approach 
is of particular importance here as we need to consider 
rather large system sizes to get accurate results for a^y 
Other methods such as the real-space formula for the 
Ghern number based on projected position operators^^AS 
(which was also used in Ref. TO) involves projection op¬ 
erators onto the occupied states, and hence requires a 
full diagonalization the 3D Hamiltonian. The diagonal- 
ization process is memory intensive and very time con¬ 
suming leading to long processing times which scale as 
0{{LxLyLz)^) compared to 0{Lj;{LyL^)^) for the Green 
function method. 

To calculate axy with high accuracy, we choose the sys¬ 
tem size by increasing the lengths in all directions until 
a^y converges (Eq. (2.3)) and further increases do not 
modify the result significantly. Our observation is that 


one usually needs to go to larger lengths (> 30) in the 
X and y directions while rather small lengths (^ 10) in 
the longitudinal z-direction often work well. This can 
be understood by noting that the system needs to be 
large enough in the xy plane to yield enough k-points for 
precise evaluation of the Ghern number, and should be 
long enough in the z-direction to give a sufficient num¬ 
ber of layers between the Weyl nodes. For calculating 
axy in this paper, the boundary conditions (BG) across 
each layer are chosen to be periodic along the z direction 
and open along y direction. We note that while axy in 
the WSM and NI phases does not depend on BG in z 
direction, it does depend on the BGs in case of 3D TIs 
(studied in Sec. H G) due to the extra half-integer surface 
Hall conductance arising from the gapped Dirac surface 
states. 


A. Stack of Ghern insulators C=1 


As the simplest model of the WSM, we consider the 
two-band Hamiltonian of multi-layer GI modeU®d7 


[clc+e,S^tas-ra3)c^ + h.c. 

X 

s=l,2 

“ f [ci+,,3(T3Cx-kh.C. 

X 


-|- m E C^O-sCx 

X 


(2.4) 


where c, denote two component fermion operators, ai 
are Pauli matrices representing spin, and are cubic 
lattice unit vectors |as| = 1. The first two sums de¬ 
scribe a stack of GI layers coupled to their neighboring 
layers through the third sum with a tunneling amplitude 
fo. Fourier transforming to momentum space, the Bloch 
Hamiltonian is given by 

/i(k) =tsuikx ai +t sin ky a2 

+ (m — r cos kx — r cos ky — t^ cos kz)a3. (2.5) 

We fix the parameters r = t = 1, and focus on the regime 
0 < ts < TO. The phase diagram is symmetric under 
TO —>■ —TO up to flipping the sign of axy In the limit 
TO —)■ 00, this model realizes an NI. As to is tuned down to 
zero, it supports two qualitatively distinct WSM phases: 
WS r for TO < ta and WS IT for |to— 2r| < t^. The WS I 
phase contains two pairs of Weyl nodes at (0, tt, ±fco) and 
(tt, 0 ,±A:o) yielding axy = l — 2ko/'K in which cos“^ m/t^. 
The WS H phase hosts a pair of Weyl nodes at (0, 0, ±fco) 
and has axy = kofn where fco = cos“^(to — 2r)/t3. The 
intermediate region between these two WSM phases is 
described by a 3D GI (i.e., a 3D T-breaking weak topo¬ 
logical insulator) where axy = I- 

Figure 1 presents the phase diagrams in the (to, W) 
plane for two different values of the interlayer hopping 
fo = 0.5 (top) and fo = 1.0 (bottom). In the former case, 
two WSM phases (denoted by WS I and WS H) are sep¬ 
arated by a gapped 3D GI phase (i.e., a 3D T-breaking 
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FIG. 1. (Color online) Phase diagram of the layered Chern 
insulator model given in Eq. (2.4). The color map represents 
Uxy and the solid lines representing phase boundaries are de¬ 
termined by the localization length. We show the phase dia¬ 
gram for two different values of ta: fa = 0.5 (top) and fa = 1 
(bottom). The system size is 30 x 30 x 12. Black triangles on 
the i-axis mark the clean phase boundaries. The dashed line 
in the lower panel connects maxima of a^y as a guide for the 
eye. We note that there is no difference in meaning between 
the solid black and solid yellow lines, the color difference is 
just to add contrast. 


weak TI phase). As the colormap shows, in both cases 
cfxy can be used to distinguish various phases even at 
finite disorder. The solid lines indicate the phase bound¬ 
aries determined by the scale invariant points where the 
localization length ^ is size independent. Comparisons of 
Gxy and ^ for few values of m are shown in Fig. 2. The top 
panels in Fig. 2 also compare Uxy for different sizes and 
confirm that the finite size corrections are rather small 
once we have reached a certain size. 

There are several interesting features in these plots. 
As the system becomes diffusive, denoted by the label 
diffusive metal (DM), axy starts to decrease and even¬ 
tually vanishes for ultra-strong disorder {W > 16) in 
the Anderson localized phase. In addition, the WSM 
phases tend to broaden as the disorder becomes stronger. 
The critical points in the clean limit separating WS I- 
3D CI-WS II-NI phases (Fig. l(top)) or WS I-WS II-NI 
phases (Fig. 1 (bottom)) are moved to the right as disor¬ 
der strength increases. In particular, axy increases when 




(b)m= 1.0 



W 



FIG. 2. The Hall conductivity and localization length ver¬ 
sus disorder strength for layered Ghern insulator model with 
fa = 0.5. In the lower panel of each subfigure, the col¬ 
ors show different cross-section sizes L = 8(blue), lO(red), 
12(orange), and 14(cyan). The dashed vertical lines repre¬ 
sent phase boundaries, the dashed horizontal lines indicate 
the value of a^y in the clean limit using the analytical expres¬ 
sions explained in the text, and the system size for a^y are 
30 X 30 X 12 (red), 70 x 70 x 16 (blue), and 90 x 90 x 20 (green). 


disorder increases in the WS II phase (Fig. 2(c)). 

This behavior, and many of the results that follow, 
can be understood by thinking about the clean limit sub¬ 
ject to weak disorder. Since our localization length data 
shows that, starting in the WSM phase, there is no crit¬ 
ical behavior until the DM phase is reached, we might 
hope that our weak disorder understanding will apply for 
a wide-range of disorder strengths; indeed by comparison 
with our numerical results the application appears suc¬ 
cessful. Now, based on the picture of a WSM as Cl layers 
in momentum space^®’^^’'*^, then as far as the low-energy 
Hamiltonian is concerned, the 3D BZ can be viewed as 
layers of Cl with a varying mass parameter 

/i(k) « t(kxai + kya2) + {M(k:,) + ^^k\)a^ (2.6) 

where = kx + ky is the in-plane momentum, M{kz) = 
m — 2r — coskz, and mj_ = r. Hence, in the direction 
of separation between the nodes the BZ can (usually) 
be divided into topological layers M{kz) < 0 of non-zero 
Chern number and trivial layers M(kz) > 0 of zero Chern 
number. The planes containing the Weyl nodes at kz = 
i/co we will call the ‘Weyl planes’ where M(±fco) = 0. 
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(a) C = 0 C = ±1 C = 0 




(c) 


y 



X 


z 


C = ±l C = =Fl C = ±l 



FIG. 3. The effective shifting directions of the Weyl nodes 
inside BZ for the layered Cl model as the disorder strength is 
increased (but below the diffusive metal regime). The Chern 
number of the momentum slices in each section of the BZ is 
shown above them. Weyl cones of different colors represent 
opposite chiralities. This motion is determined by the SCBA 
calculation and is consistent with the numerical results. For 
subfigure (c) the Weyl nodes do not move. 


With this set up, let us consider the effects of disor¬ 
der. Consider two Weyl planes in the BZ where either 
the layers inside (Fig. 3(a)) or outside (Fig. 3(b)) the two 
planes are topological and the complementary region is 
trivial. Numerically, our results can be consistently in¬ 
terpreted if the effect of disorder moves the Weyl planes 
in the direction such that the number of topological lay¬ 
ers increases (Fig. 3). This means that insulating layers 
near the Weyl planes on the trivial side become topologi¬ 
cal with non-zero Chern number. This effect can happen 
since the disorder can be thought to renormalize the mass 
parameter down to negative values and hence inverting 
the bands. This type of behavior has also been shown for 
a low-energy model Cl using the SCBA and including the 
quadratic (in momentum) mass terms^^. Using only the 
Born approximation (i.e., not self-consistent), they found 
the mass renormalization 5M to be^"^ 


5M = - 


487rm_L 


In 


7r^771_L 

2M 


(2.7) 


One important consequence of this result is that sign of 


SM depends on the sign of quadratic term m±k'j_ as fol¬ 
lows: if m_L > 0, then 6M < 0 and vice versa. 


For our model, m± = r is always positive and hence 
the mass renormalization is always negative. This is then 
a way to understand the underlying process leading to an 
enhanced AHE response in the WS II phase as a function 
of disorder strength. Moreover, we observe that the insu¬ 
lating phase NI in the vicinity of WS II transitions into 
a WSM at a finite disorder before transforming to the 
DM (Fig. 2(c)). This can also be attributed to the nega¬ 
tive mass renormalization which makes the trivial layers 
near the origin topological and lets Weyl nodes nucleate 
at the interface between newly formed topological layers 
and trivial ones. 


However, a different situation occurs in the WS I 
phase; i.e., axy starts as a plateau with initially increas¬ 
ing disorder, and then decreases as disorder is further 
increased (Fig. 2(a)). In this phase, there are two pairs 
of Weyl nodes with the same chirality in the two planes 
at kz = ±fco allowing for the Chern number to change 
by two as the Weyl planes are crossed in momentum 
space (Fig. 3(c)). Hence, the layers between Weyl planes 
\kz\ < fco have C = —1 and the others \kz\ > k^ have 
(7 = 1. Hence, all of the layers are topological in the 
clean limit. This is different than the usual form seen in 
the WS H phase. 


Interestingly, unlike WS H, we find that the SCBA 
mass renormalization for these nodes is zero, and hence, 
the trend of in this case cannot be immediately deter¬ 
mined as a perturbative effect as it can be in WS H. This 
is one reason why we might expect (and as we do see nu¬ 
merically) that axy should remain constant at weak disor¬ 
der. To help understand this behavior beyond perturba¬ 
tion theory, we numerically calculated the phase diagram 
of the 2D Cl in Appendix A. By applying the results of 
this phase diagram we expect that rather weak disorder 
can localize the edge modes of Cl layers on both sides 
of Weyl planes, which have opposite Chern number, and 
will nucleate a region with (7 = 0 beginning in each Weyl 
plane. One might understand this by considering that, 
since nearby momentum slices on the opposite sides of the 
Weyl plane have opposite Chern number, they can couple 
at weak disorder and annihilate. In fact, since, the reduc¬ 
tion of Chern number on both sides would be symmetric, 
we would again expect that we will not find any change 
in axy initially. However, we note that fco < '^/2 which 
means the region with (7 = 1 {\kz\ > fco) is wider than 
the region with (7 = —1 {\kz\ < fcg ). Thus, as we con¬ 
tinue increasing the disorder, the central regionj/c^l < fcp 
becomes entirely trivial first, while part of the \kz\ > kg 
region remains topological. So, a further increase of dis¬ 
order results in making these layers trivial and conse¬ 
quently reduces axy 
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B. Stack of Chern insulators C=2 

A simple model of a double Weyl semimetal (dWSM) 
can be constructed by stacking layers of CIs with Chern 
number two, 


^ =\Y1 [4+a,(io'i - 


Cx + h.c. 


m 




s=l,2 

^ 2 ^ \ 

X 

s=l,2 


X [4+aV^2Cx - 4+^,^ azCx + h.C. 


2 


--Eh 


x+a 3 


CTaCx + h.C. 


( 2 . 8 ) 


where are unit vectors for the cubic lattice, and a(j_ = 
ai ± a 2 connect next nearest neighbors in the xy-plane. 
In the clean limit, the Bloch Hamiltonian can be written 


as 


h{k) =t{coskx — cosky) cti + 
+ (m 


. 3 , CUsAiy) fX 

r cos kx — r cos ky 


2t' sin kx sin ky a-i 

sky — cos kz)<Js- 


(2.9) 


When |to — 2r| < this model yields Weyl nodes at 
(0,0,±fco) each with a Berry-monopole charge of ±2 
where fcp = cos“^(to — 2 r)/t3. The Hall conductivity 
is given by Oxy = 2ko/TT. 

The phase diagram for = 2r is shown in Fig. 4. It 
essentially follows the same properties of the multi-layer 
Cl model in the previous section with (7 = 1. Similar to 
that model, in the diffusive limit, the Hall conductivity 
decreases monotonically to zero. Moreover, before reach¬ 
ing the diffusive limit we can understand the structure of 
the phase diagram using the Born-approximation argu¬ 
ments from the previous section. The clean Hamiltonian 
for a layer at and near the origin of the in-plane mo¬ 
menta kj^ = (kx, ky), reads 


h{k) - {kl- kl)ai + 2t'kxkya2 


t 

2 ' 

{M{kz) -I—^fci)o’3, 


( 2 . 10 ) 


in which M{kz) = m — 2r — t^coskz determines trivial 
layers M(A; 2 ) > 0 from topological ones ^/(A:^) < 0. Since 
the phase boundary of the dWSM moves further to the 
right (larger m) then we need to understand the instabil¬ 
ity of trivial layers near the Weyl nodes, M(±fco) = 0, to 
transition into the topological phase with Chern number 
two. Nicely, this can again be explained in terms of mass 
renormalization towards the topological phase similar to 
the (7=1 case^^. For t = 2t', the SCBA, in the Born 
approximation limit, can be analytically derived to be 


SM = - 




m± 


+ t^) 

4M2 


96tt 




In 


2 M 


( 2 . 11 ) 


S 4 



This implies that the Weyl nodes move away from one an- 


FIG. 4. (Color online) Phase diagram of the dWSM model 
given by Eq. (2.8). The color map represents a^y, and the 
solid lines are phase boundaries determined by the localization 
length. The black triangle on the m-axis marks the clean 
phase boundary. The system size is 30 x 30 x 14. We note that 
there is no difference in meaning between the solid black and 
solid yellow lines, the color difference is just to add contrast. 


other towards the BZ boundary so that axy will increase 
(similar to Fig. 3(a)). 

Let us now comment on the similarities between the 
WSMs realized by multi-layer Cl models. In terms of 
the Hall conductivity, they share two common proper¬ 
ties: (i) in the weak disorder limit the Weyl nodes are 
pushed toward the region with zero Chern numbers (see 
Fig. 3), hence expanding the topological portion of the 
BZ and (ii) in the diffusive limit the AHE is weakened as 
the disorder is made stronger. It is worth noting that the 
former property is specihc to the multi-layer Cl systems 
where the sign of quadratic term rk'^ is independent of 
the position ko of the Weyl nodes. The model in the next 
section that is constructed from a 3D Dirac model, will 
not always show this behavior in the entire phase diagram 
since the sign of quadratic term depends on the nodal 
momenta. Given our localization length data, we have 
not noticed any drastic differences between the dWSM 
and the ordinary WSM in terms of their stability against 
disorder or their transitioning to the DM phase at fi¬ 
nite disorder. This is in contrast with the conclusions of 
Ref. 69, though we are not sure if the model used in that 
work should be directly comparable to ours. If they can¬ 
not be directly compared then perhaps this would resolve 
the discrepancy. 

In terms of the longitudinal conductivity, our calcula¬ 
tion of localization length suggests that the WSM stays 
semi-metallic at weak disorder and subsequently becomes 
a DM for sufficiently strong disorder. In other words, it 
never becomes insulating in this process. As we will see 
in the next section, the WSM obtained by breaking the 
time-reversal symmetry in the 3D Dirac system does not 
follow this trend and transitions to an intermediate insu¬ 
lating phase along the way transitioning into a DM. 
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C. 3D Dirac semimetal subject to a Zeeman field 


Consider the four-band model of a time-reversal invari¬ 
ant Z 2 TfSO’Si 


H =- 4-Sa,(*tos-»’/3)Cx+ll.C. 


s=l,2,3 

I- (m -I- Sr) ^2 4/3cx, 


( 2 . 12 ) 


where the Dirac matrices are given by 


Cks = D C) (Ts 


0 CTs \ 

a, 0 J’ 


/3 = T 3 (g) 1 = 


1 0 \ 

0 -1 )' 


In this convention the cr and r matrices act on the spin 
and orbital degrees of freedom respectively. We set the 
parameters at r = t = 1 from now on. 

The Bloch Hamiltonian in the clean limit is: 


^(k) = 

s=l,2,3 


tag sin kg 


rj3 cos kg 


(to -I- 3r)/3. 


There are two important symmetries of this model: 
(i) time-reversal symmetry (TRS) with the operator 
T = icF 2 lC such that tT 2 /i(k)o ’2 = /i*(—k); (ii) inver¬ 
sion symmetry (IS), represented by X = T 3 X, such that 
T'ih{k)T'i = h{—k). This model realizes a strong TI (STI) 
phase over the mass range — 2 r < to < 0 , and there is a 
transition from STI to an NI at to = 0. The critical point 
at TO = 0 is described by a gapless 3D Dirac Hamiltonian 
with four bands touching at a single point. If we add a 
TRS or IS breaking term then the critical point broad¬ 
ens into an intermediate semi-metallic region, i.e. a Weyl 
semimetal. In this case the fermions with opposite hand¬ 
edness are shifted away from the T-point and form Weyl 
nodes at different points inside the BZ. 

One such term we could add to break TRS is a Zeeman 
splitting field: 


Hb, =bg'22 4'^oo-sCx (2-13) 

X 



m 


FIG. 5. (Color online) The phase diagram of a WSM gener¬ 
ated by broken TRS in a 3D Dirac model given by Eq. (2.12). 
The color map represents aj;y and solid lines are phase bound¬ 
aries determined by localization length. The AHI phase 
is contained in a closed region of the phase diagram and 
the NI-AHI and AHI-DM boundaries eventually intersect at 
(m, W) = (2.6, 9.5) (this lies outside the graph on the top 
right). The black triangles on the m-axis mark the clean 
phase boundaries. 63 = 0.6 and the system size is 30 x 30 x 10. 
We note that the STI phase would have a quantized surface 
Hall conductivity if we had chosen open boundaries in the 
^-direction instead of periodic. We checked that this was the 
case, but did not include both diagrams. 


towards the BZ boundary until TOc = — t < 

0 where they reach their maximum distance. Making 
TO smaller results in bringing the nodes closer and they 
finally annihilate each other at the origin once to = —63 
is reached. After annihilation the STI phase is formed. 

The spectrum near a Weyl point at (0,0, fco) is given 
by 


h{k) = v±{kxrii + kvV2) + {vzkz + ^fci)??3, (2-15) 

where r]i are Pauli matrices in the projected space of the 
two touching energy bands, the in-plane momentum is 
= kl + ky, and the other coefficients are 


where 0 < bg < t. Depending on s = 1,2,3 this term 
splits the Dirac node into two Weyl points which are 
shifted from the origin in opposite directions (±fco) along 
the s-axis in momentum space. The amount of shift is 
given by 


cosfco = 


+ {m + t)^ — 
2t{m +1) 


(2.14) 


Therefore, in the presence of this term, a WSM phase 
forms in the range \m\ < bg between the STI and NI 
phases. From now on, let the Zeeman field be in the z- 
direction with 63 ^ 0. From the NI side (to > 63 ), as 
TO is decreased the two Weyl nodes start to nucleate at 
the origin when to = 63 and move apart from each other 


m + t — t cos ko 

V± = TOj_ = -;- 1 

h 

m + t . 

Vz = —^-csin^o. 

03 


(2.16) 

(2.17) 


Note that the chiralities of Weyl nodes at the two points 
( 0 , 0 , ±fco) are opposite to each other since only the third 
component of velocity flips sign v^i—kQ) = —v^iko). 

Figure 5 presents the phase diagram of this model 
where a^y is shown as the background color. The solid 
lines are phase boundaries determined by the scale invari¬ 
ant points of the localization length. Typical results for 
the Hall conductivity and localization length as a func¬ 
tion of disorder strength for four to points are also shown 
in Fig. 6 . For comparison, the phase diagram of Z 2 STI 

















(a) m= - 0.9 


(a)m= 0.1 



0 5 10 

W 



W 




W W 


FIG. 6 . The Hall conductivity and localization length versus 
disorder strength for the WSM generated by broken TRS in 
a 3D Dirac model with 63 = 0.6. The dashed vertical lines 
represent phase boundaries and the dashed horizontal lines 
indicate the value of in the clean limit using Eq. (2.14). 
For a^y, the system size is 30 x 30 x 10 (red) and 50 x 50 x 16 
(blue). For ^/L, the colors show different cross-section sizes 
L = 6 (blue), 8 (red), lO(orange), and 12(cyan). We note that 
in subfigure (a) the STI phase would have a quantized surface 
Hall conductivity if we had chosen open boundaries in the 2- 
direction instead of periodic. We checked that this was the 
case, but did not include both diagrams. Also, CR in panel 
(b) and (c) refers to the critical STI phase near the critical 
line between STI and AHI phases, where the non-zero bulk 
a^y could be due to finite-size effects (i.e. cr^y varies smoothly 
from AHI to STI). 


in the absence of the Zeeman field has been studied pre- 
viously®®’®^^®"^. The apparent shift of the phase bound¬ 
aries to the right (toward larger m) is controlled by the 
higher order terms in momenta in the effective theory 
Eq. (2.15). This has also been reported in the previ¬ 
ous studies and leads to the so-called disorder-induced 
3D TI phase®^’®®. However, the phase boundary between 
the STI and NI, which was originally a line®^’®®, now 
broadens into a semi-metallic phase. Notice that the top 
left region of the phase diagram in Fig. 5 with negative 
axy corresponds to the tail of another WSM phase which 
originally forms in the mass range |m -I- 2r| < 63 between 
the STI and a weak TI in the clean limit and is shifted to 
the right in the diffusive limit (see also Fig. 6 (a)), similar 
to the WSM (|m| < 63 ) that is our focus here. 

Let us remark that the label STI may be a bit impre- 



W 



4 5 6 7 8 9 10 

W 

FIG. 7. The localization length as a function of disorder 
strength in the WSM generated by TRS-broken 3D Dirac 
model with &3 = 0.6 for m = 0.1 (a) and m = 1.8 (b). 
The dashed vertical lines represent phase boundaries. The 
colors show different cross-section sizes L — 6 (blue), 8 (red), 
lO(orange), and 12(cyan). In panel (a), we should note that 
there is a hint of the scale-invariant point between AHI-STI. 


cise in presence of TRS breaking Zeeman term. However, 
in what we call STI, the Zeeman term 63 is weak com¬ 
pared to the bulk gap m so that the bulk wave-functions 
are not affected much and can be adiabatically connected 
to those of the time-reversal symmetric STI in absence 
of the Zeeman term. This implies that the bulk STI (in 
presence of weak Zeeman) has a zero Hall conductance. 
However, the surface states are gapped by the Zeeman 
term and this leads to surface AHE. In Fig. 6 (a), we im¬ 
pose periodic BCs along the z-direction and the surface 
contribution is absent hence leading to a^y = 0. For 
consistency we checked that open BCs in the z-direction 
yields the quantized AHE {a^y ■ = 1) due to the 

gapped Dirac surface states, although we do not present 
the figure here. The latter set-up is very similar to a 
time-reversal symmetric STI along with TRS breaking 
ferromagnetic layers on the top/bottom surfaces. Conse¬ 
quently, the STI phase is clearly different from NI, where 
there is neither bulk nor surface Hall conductance. In the 
clean limit it can also be thought of as a TI protected by 
inversion symmetry. The key feature is that the system 
still exhibits the characteristic topological magnetoelec¬ 
tric effect (quantized surface Hall conductance). 
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In addition, we surprisingly find that the semi-metallic 
phase is interrupted by an insulating phase before turning 
into the DM. This insulating phase -which we call a dis¬ 
order induced anomalous Hall insulator (AHI)- must be 
contrasted from the NI and STI as it has a non-zero bulk 
Hall conductivity (the NI does not have a Hall conduc¬ 
tivity and the STI, at most, has a surface Hall conductiv¬ 
ity). The AHI emerges in the strong disorder regime as 
a result of localizing the Weyl nodes without annihilat¬ 
ing each other; i.e., in the effective BZ picture, the nodes 
are far away from each other and localize on their own. 
Therefore, the AHI phase with non-quantized 3D Hall 
conductivity (per layer) axy (as opposed to the weak- 
disorder limit of the 3D Cl studied above where a^y = 1 
per layer) is an exotic insulating disordered phase which 
cannot be realized in systems with translational symme¬ 
try (which implies the a^y is quantized to be the same 
multiple of e^/h for each layer®®’®®). It is important to 
note that from our numerics that if we view the AHI as 
a quasi-2D insulator, it could be interpreted as a Hall 
liquid with a'^y ^ 0, but there is another possibility 
that the AHI phase in the thermodynamic limit con¬ 
verges to the “Hall insulator” state, originally introduced 
in Ref. 87, and also experimentally observed recently®® 
at the magnetic-field-driven superconductor to insulator 
transition of amorphous indium oxide, where —>■ 0 

and a^y —>■ 0, in such a way that a'^y oc Y so 
that Pxy becomes finite. In order to fully characterize 
the AHI phase one needs to compute a^x and a^y in the 
thermodynamic limit at low temperatures and we cannot 
rule out either possibility based on our current numerical 
data. 

Based on our data, the physical difference between STI 
and AHI is the existence of chiral edge modes (in the xy 
plane circulating around the z axis) in the latter when 
the boundary condition is periodic along the z direction. 
In other words, the chiral edge modes are descendants 
of the chiral Fermi arcs in the Weyl semi-metallic phase 
that the AHI is originated from. We anticipate that in 
order to transition from the AHI phase to either of STI 
and NI phases, the system must undergo a (mobility) 
gap closing (scale-invariant critical point) where the ex¬ 
tended bulk states emerge, leading to a sharp change 
in the topological response (Hall conductivity) in these 
insulating phases. Our localization length analysis sug¬ 
gests such a behavior at the NI-AHI and STI-AHI phase 
boundaries (see Fig. 7). As mentioned in the caption of 
Fig. 5, we checked that the AHI phase occupies a closed 
region of the phase diagram by finding the intersection 
between the converging phase boundaries of the DM-AHI 
transition and the NI-AHI transition. 

The last remark we want to make is about the di¬ 
rection in which the Weyl nodes are shifted by weak 
disorder. Our finding based on axy implies that for 
m > rric (defined above) the Weyl nodes tend to move 
away from one another leading to an enhancement of 
the AHE (Fig. 6(c), (d)); while for m < rric the Weyl 
nodes are moved closer together so that they annihilate 


at finite disorder (Fig. 6(b)). In both cases, the layers 
with non-zero Chern number are located between two 
Weyl nodes, however, the nodes can move either towards 
or further from each other. This is different from the 
multi-layer Cl realizations in previous sections where dis¬ 
order always pushes the nodes apart to create a wider 
topological region (Fig. 3(a)). We can understand this 
for 3D Dirac model by realizing that the coefficient of 
the quadratic term to_l in Eq. (2.16) changes sign as fcp 
crosses the critical value kc = sm~^{b^/t) corresponding 
to the critical mass rric- This causes the mass renormal¬ 
ization (Eq. (2.7)) to pick up the opposite sign depending 
on whether m is larger or smaller than rric'. in the for¬ 
mer, nodes are moved apart (similar to previous sections) 
whereas in the latter, nodes are pushed towards one an¬ 
other. Thus, the movement direction of Weyl nodes in 
this case coincides as well with SCBA analysis as long as 
this complication is taken into account^^. 


III. DISCUSSION 

In conclusion, we have studied three different lattice 
model realizations of Weyl semimetals in the presence of 
an on-site disorder potential by calculating the localiza¬ 
tion length and the Hall conductivity. The first two mod¬ 
els are based on stacking layers of a 2D Chern insulator 
with Chern numbers one and two respectively, whereas 
the third model is obtained by breaking the TRS in a 
3D four band Dirac model. We have found that in a 
WSM phase, weak disorder may effectively act to cause 
the Weyl nodes to move inside the BZ. The movement 
direction is not determined by the position of the Weyl 
nodes, e.g., whether they are close to the origin or bound¬ 
aries of the BZ (as opposed to Ref. 70). Instead their 
movement is actually determined by the Chern number 
of the 2D (momentum-space) layers in the clean limit 
where the layered structure of the 3D BZ can be applied. 
In short, disorder forces the Weyl nodes to expand into 
the region with zero Chern number instead of the re¬ 
gion with non-zero Chern number. We have successfully 
applied this idea to interpret all our results. Moreover, 
we have also shown that the SCBA analysis is consistent 
with this picture in the case of WSMs with two Weyl 
nodes. We have seen however, in the case of a WSM with 
four nodes where all layers carry non-zero Chern num¬ 
ber, the SCBA formula predicts no motion for the Weyl 
nodes although our data clearly shows a different phe¬ 
nomenon. Remarkably, this effect can be explained with 
the aid of the phase diagram of a single disordered 2D 
Chern insulator, which implies that layers with opposite 
Chern numbers in the vicinity of Weyl nodes hybridize 
and start to lose their Chern number. Nevertheless, an 
explicit perturbative understanding of this trend remains 
an open issue. 

We have also observed that weak disorder may nucleate 
Weyl nodes within the BZ of insulating phases near the 
WSM phases hence endowing them a finite Hall conduc- 
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tivity; or oppositely, weak disorder may annihilate Weyl 
nodes in the WSM phases near insulating phases leading 
to a TRS broken insulator adiabatically connected to a 
weak or strong topological insulator. We have found that 
both types of transitions are consistent with the predic¬ 
tions of SCBA. 

In addition, we should note that our AHE response 
data in the WSM phase is consistent with the recent an¬ 
alytical field-theoretic results^^’®®, once the underlying 
assumptions of their calculations are taken into account. 
These studies have shown that the AHE response of a 
WSM is robust against weak disorder and always given 
by the separation of the Weyl nodes (which is an intrinsic 
property of the WSM). As a result, they found that axy 
remains intact even in the diffusive limit. These results 
rely on two important assumptions: first, the higher or¬ 
der (quadratic at lowest) corrections to the band disper¬ 
sion are neglected and second, the non-perturbative local¬ 
ization effects which ultimately lead to the Anderson lo¬ 
calization are not included. The first assumption implies 
that the distance between Weyl nodes is not renormalized 
by the disorder and the second assumption prevents any 
process causing axy to decrease. However, these two as¬ 
sumptions cannot be made in our numerical calculations 
and hence the apparent difference between our numerical 
analysis and the continuum field-theoretic calculations 
are unavoidable. Nevertheless, as we show here, the ob¬ 
served trend of AHE in the WSM can be related to the 
renormalized distance (rather than the original distance 
in the clean-limit) between two Weyl nodes and the fact 
that AHE in the diffusive limit decays to zero slowly is 
consistent with the coexistence of diffusive transport and 
AHE as argued in the continuum limit calculations. 

Finally, our localization length data signals a surpris¬ 
ing effect that is special to the 3D Dirac model and was 
not found in the multi-layer systems: the WSM region of 
the phase diagram is separated from the diffusive metal 
through an emergent Hall insulator with a non-quantized 
(per layer) bulk Hall conductivity. Although the exis¬ 
tence of a disordered anomalous Hall insulator with a 
non-quantized bulk Hall conductivity is not violated by 
any physical principles, it is not really clear to us that 
why this phase only appears in the 3D Dirac model. One 
interesting direction is to possibly characterize this phase 
further and check if it survives in the thermodynamic 
limit. 
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FIG. 8. Phase diagram of a disordered 2D Chern insulator 
that serves as a building block of the 3D WSM phase with 
(a) isotropic (ri = r 2 = t) and (b) anisotropic (ri = 1.3t, 
r 2 = 0.7t) Wilson-Dirac mass terms. The chemical potential 
is set to zero. 


Appendix A: Disordered Chern Insulator 

In this appendix, we use the localization length and 
Chern number to numerically determine the phase dia¬ 
gram of a disordered 2D Chern insulator model (Fig. 8 ). 
We use a generic form of the two band model, 

h(k) =tsiTikx ai +t sin ky a 2 

+ {m — ri cos kx — r 2 cos ky)az- 

where the Wilson-Dirac mass terms cos ki can have 
anisotropic coefficients. The WSM studied in Sec. HA 
consists of isotropic layers of this model where ri = r 2 = 
r. As Fig. 8 (a) shows, the clean limit has critical points 
at TO = —2r and m = 2 r that are bent towards the NI 
phases with increasing disorder. This is consistent with 
the SCBA formula (Eq.(2.7)). Interestingly, the clean 
critical point at to = 0, separating the two Cl phases, 
splits to form two phase boundaries and a trivial insu¬ 
lating phase appears in-between them. This behavior is 
beyond the perturbative SCBA, which gives an exactly 
vanishing mass correction in this region. In fact, the 
SCBA integrals for the continuum limit model are van¬ 
ishing in this case due to isotropy of the system. From 
a low-energy perspective, at to = 0 the isotropic system 
carries simultaneous 2D Dirac nodes at (0, tt) and (tt, 0) 
in the clean limit. Then, at finite disorder, each of the 
new phase boundaries must essentially contain one Dirac 
node to allow for a change in the Chern number by one. 

In other words, as we tune to from -|-oo to —oo at a 
finite disorder (e.g., at lU > 2 ), we should encounter four 
critical points to go through the NI-CI-NI-CI-NI sequence 
of phases. Such a process is possible at finite disorder, 
since the two band simultaneous crossings to = 0 are 
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protected by the C 4 lattice symmetry which is lifted in 
the presence of disorder. For each disorder realization 
this symmetry is explicitly broken and in principle, one 
crossing may occur earlier than the other. However, on 
average, the C 4 symmetry is recovered and the crossing 
at each boundary cannot be attributed to only one node 
at ( 0 , 7 r) or ( 7 r, 0 ). In fact, we checked that this is the 
case, by looking at the Fourier transform of the eigen¬ 
states (near zero energy) of the disordered Hamiltonian 
on either phase boundaries. We found that the averaged 
peak heights at (0 , tt) and (tt, 0 ) are equal. 

Note that the above situation, where the SCBA van¬ 
ishes, is specific to the Ci symmetric model. It is easy 
to see that if one explicitly breaks the C 4 symmetry in 
the clean limit by choosing ri different from r 2 , the two 
band crossings do not occur simultaneously and there 
will be an intermediate trivial insulator between two CIs 


for |m| < |ri — r 2 | (Fig. 8 (b)). Remarkably, the SCBA 
integral in this case would be non-zero, proportional to 
the anisotropy dm oc ±|ri — r 2 |/(ri -|- r 2 )^, and have a 
correct sign consistent with the curvature of the phase 
boundaries given by the numerics. Moreover, we ob¬ 
served that the momentum distribution of the energy 
eigenstates (near zero-energy) on the phase boundaries 
in the strongly disordered limit (IF > 6 ) of both the 
isotropic and anisotropic models have equal weights at 
(0, tt) and (tt, 0), meaning that the original difference in 
the clean limit between these two models is completely 
lost as a consequence of large scattering between various 
momentum modes. This is opposite for the phase bound¬ 
aries at weaker disorder where the anisotropic model has 
momentum weight dominated by (0, tt) or (tt, 0) whereas 
in the isotropic case, as mentioned above, the averaged 
peak heights of both momenta are the same. 
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